A full Eulerian finite difference approach for solving fluid-structure coupling problems 
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A new simulation method for solving fluid-structure coupling problems has been developed. All 
the basic equations are numerically solved on a fixed Cartesian grid using a finite difference scheme. 
A volume-of-fluid formulation (Hirt & Nichols (f98f, J. Comput. Phys., 39, 201)), which has been 
, widely used for multiphase flow simulations, is applied to describing the multi-component geometry. 

The temporal change in the solid deformation is described in the Eulerian frame by updating a 
left Cauchy-Green deformation tensor, which is used to express constitutive equations for nonlinear 
Mooney-Rivlin materials. In this paper, various verifications and validations of the present full 

■ Eulerian method, which solves the fluid and solid motions on a fixed grid, are demonstrated, and 
' the numerical accuracy involved in the fluid-structure coupling problems is examined. 

^ , PACS numbers: 47.11. Be, 47.10.A-, 47.55.-t 

^— 1 . 

I. INTRODUCTION 

PK 

f*) |' Fluid-Structure Interaction (FSI) phenomena appear in many places, e.g., biological systems, and industrial pro- 
cesses. In life science, the numerical prediction of the FSI problems involves twofold significance: fundamental and 
q ■ practical aspects. For example, it would allow us to gain insight into how life is sustained, and support surgical 
' planning in medical treatments. Conventionally, the computational fluid dynamics is modelled in an Eulerian way, 
^ . while the computational structure dynamics is normally treated in a Lagrangian way. The coupling of the fluid and 
O ' structure dynamics is a formidable task due to such a difference in the numerical framework as well as its multi-physics 

■ nature. 

^-y In a creeping (zero Reynolds number limit) flow involving a biological membrane, once boundary conditions on the 
\ membrane surface are determined via constitutive laws for e.g. an in-plane stress and a bending moment, the bulk 
flow field is found using Green's function. In such a situation, no volumetric mesh is needed in the bulk region. A 
boundary element method is applicable to predicting capsule [59] and red blood cell [60] motions interacting with the 
t— I ■ fluid flow. 

For non-zero Reynolds number flows, on the other hand, it is necessary to set out the computational mesh over the 
, entire domain. There are currently several major approaches classified with respect to the computational treatment 
' how the kinematic and dynamic interactions are coupled on the moving interface. 

The most accurate approach is raised as Arbitrary Lagrangian Eulerian (ALE) [3, 22, 29, 30, 50] or Deforming- 
• ■ Spatial-Domain/Spacc-Timc (DSD/ST) [31, 75, 76] technique, in which the body-fitted mesh is used. The latter 
method enables one to arbitrarily choose spatiotemporal nodes, and facilitates to simulate the moving particle and 
FSI problems. These approaches are referred to as an interface-tracking approach, in which the surface mesh is 
accommodated to be shared between both the fluid and solid phases, and thus to automatically satisfy the kinematic 
condition. As the flow is resolved along the moving/deforming object surface, the boundary layer is expected to be 
highly resolved, and the dynamic interaction is accurately coupled through iterative procedures. Once the body- 
fitted mesh is provided, the state-of-the-art of the interface-tracking approach is satisfactory for achieving accurate 
predictions, including the applications of moving rigid particles [25, 37], moving hyperelastic particles [12], parachute 
aerodynamics [68], blood flows [2, 11, 74, 77-79] and heart flows [84, 95]. However, the whole computational domain 
has to be re-meshed as the object moves/deforms, which is computationally expensive. 

An alternative to the interface-tracking approach is an Eulerian-Lagrangian approach, in which the fluid and solid 
phases are separately formulated on the fixed Eulerian and Lagrangian grids, respectively. A noticeable contribution 
is the development of the Immersed Boundary (IB) method by Peskin [57, 58], who introduces a pseudo delta function 
for communication between the Eulerian and Lagrangian quantities, and demonstrated the landmark simulation of the 
blood flow around heart valves [57, 58]. The Fictious Domain (FD) method [14, 15, 65, 90], and PHYSALIS [71, 96] 
for specific multiphase flow problems with circular or spherical particle are also classified into the Eulerian-Lagrangian 
approach. The IB and FD methods have been applied to a variety of studies, for example, moving rigid particles 
[14, 15, 39, 40, 91], moving flexible bodies [27, 49, 97], red blood cell motion [10, 16, 49], and restricted diffusion 
with permeable interfaces [26] . The IB method has also inspired many researchers to propose a number of improved 
methods. For example, to facilitate the application to medical problems, the immersed finite element method, in which 
a new kernel function instead of the pseudo delta function is introduced for determining the cut-off region around the 



interface on a non- uniform mesh system, is proposed in [44, 92-94]. Also, the immersed interface treatment [41-43, 73] 
improves the sharpness of the interface by incorporating the jump in the stress and velocities across the interface. 
Recently, a conservative momentum exchange method [72] is proposed to facilitate the simulation of the fluid flow 
involving a number of elastic particles by combining the finite difference and finite clement approaches. 

In the above-mentioned methods for predicting the motion/deformation of hyperelastic material, the solid displace- 
ment is temporally updated in the Lagrangian way. In general, the hyperelastic constitutive law is expressed as a 
function of the deformation gradient tensor F = dx/dX, here x denotes the current coordinates, and X the reference 
coordinates [4] . The label of each Lagrangian node links the reference configuration and the set of the current node 
positions representing the current configuration. Therefore, the Lagrangian description has been preferably employed. 
However, the re-meshing procedure at each time step leads to intensive computation if system involves complicated 
geometry of solid and/or a large number of objects [69]. 

Let us consider patient-specific simulations in a medical field. The multi-component geometry of the inner side of 
the human body is provided as voxel data, which are converted from the medical image of a Computed Tomography 
(CT) or a Magnetic Resonance Imaging (MRI). The basic idea of the medical image-based simulation is found in 
[74, 77], in which the voxel data are converted into the finite element mesh before starting the computation, and the 
mesh is reconstructed with time advancement. As pointed out in [47, 51, 88], since technical knowledge and expertise 
are required for the mesh generation and reconstruction, an alternative simulation method without mesh generation 
procedure would be desirable to extend the FSI simulation to certain additional classes of problems in the medical 
field. Motivated by such a practical background, the full Eulerian finite difference methods, which directly access the 
voxel data to describe the rigid boundary on the fixed Cartesian mesh and avoid difficulty in mesh generation and 
reconstruction, have been developed in [47, 88]. The application includes the prediction of a coil embolization for 
aneurysms [51], which has been used in practice to support surgical planning. 

We further develop a full Eulerian method for fluid-structure interactions involving flexible hyperelastic material. 
In the interface-capturing methods for multiphase flow simulations, for instance, Volume-Of-Fluid (VOF) [23], level 
set [6, 55, 56, 64, 70], and phase field [35, 85] approaches, one set of governing equations for the whole flow field, 
referred to as a one-fluid formulation [81], is often employed. We follow such an idea, and write all the basic equations 
on a fixed Cartesian grid in a finite difference form. Several Eulerian solvers for modelling the solid deformation 
have been proposed for linear elastic materials [47, 87] , for elasto-plastic materials [54, 82] for nco-Hookean materials 
[9, 46, 83] and for more general hyperelastic materials [8]. We treat a Mooney-Rivlin hyperelastic material [48, 61], 
and consider the nonlinearity of the Cauchy stress with respect to a left Cauchy-Green deformation tensor [4, 20]. 
Our fluid-structure coupling approach is characterized by the feasibility in implementing the hyperelastic constitutive 
law into the standard incompressible fluid flow solvers, and also by the treatments in capturing the fluid-structure 
interface and the solid deformation. In the Lagrangian method, the body-fit finite elements automatically distinguish 
between the fluid and solid phases. In the present Eulerian method, considering that the voxel data contain the 
volume fraction of fluid and solid, we apply the VOF formulation [23] to describing the multi-component geometry 
(see figure 1). The large deformation is usually described by using the Piola-Kirchhoff stress tensor as a function of the 
deformation gradient, which is suited to the Lagrangian formulation as mentioned above. By contrast, the Eulerian 
formulation lacks of the material points to link between the reference and current configurations. Therefore, we must 
devise a method to quantify the level of deformation. To this end, we introduce the left Cauchy-Green deformation 
tensor B(~ F ■ F T ) defined on each grid point, and temporally update it (see figure 2). We will devote several test 
computations to the verification and validation, and investigate the numerical accuracy involved in the fluid-structure 
coupling. 

The paper is organized as follows. In §11, we present the details of the basic equations of the system consisting 
of Newtonian fluid and hyperelastic material, and formulate constitutive equations suited to the full Eulerian FSI 
simulation. In §111, we explain the simulation methods in terms of the time-stepping algorithm and the finite differ- 
ence descriptions. In §IV, to explore the validity of the advocated numerical procedure, we perform three kinds of 
tests. Firstly, we address a one-dimensional problem of the oscillatory parallel fluid-solid layers. Secondly, we make 
comparisons with available simulations. Thirdly, we examine how the shape reversibility of the hyperelastic materials 
is reproduced. In §V, we provide some perspectives to conclude the paper. 

II. BASIC EQUATIONS 

A. Governing equations and fluid-structure mixture representations 

Figure 3 shows the notation of the fluid-structure systems to be addressed. Let us consider an incompressible 
hyperelastic domain Sl s submerged in an incompressible fluid domain ttf, which is bounded with rigid flat walls. 
Hereafter, the suffices / and s stand for the fluid and solid phases, respectively. We focus on the system, where the 
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FIG. 1: Schematic figure explaining the difference between the interface recognitions of the Lagrangian and Eulerian approaches. 
In the Lagrangian method, the body-fit mesh distinguishes between fluid and solid phases, while in the present Eulerian method, 
the solid volume fraction tf>s does. The contour at <f) s — 1/2 indicates the interface. 
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(a) Lagrangian approach 



(b) Present Eulerian approach 



FIG. 2: Schematic figure explaining the difference between the solid deformation descriptions of the Lagrangian and Eulerian 
approaches. In the Lagrangian method, the relative displacement between adjacent material points from the reference to 
current configurations quantifies the deformation level. In the present Eulerian method, to quantify the deformation, the left 
Cauchy-Green deformation tensor B is introduced in the Eulerian frame, and temporally updated. 



walls are in contact with only fluid at the boundary Tyy, and the moving wall drives the fluid and solid motions. Both 
fluid and solid are homogeneous, i.e., the material properties are uniform inside each phase. We shall restrict our 
attention to the kinematic and dynamic interactions at the fluid-structure interface T/. The fluid and solid densities 
are assumed to be identical (pf = p s = p) as in many analyses for biological systems (e.g. [79, 97]), and no external 
body force (b = 0) is to be exerted on the continua. The extensions to the non-identical density and the non-zero 
body force would be readily achieved. 

For incompressible fluid and solid, the governing equations consist of the mass and momentum conservations: 

V-v f = 0, xefif, 
v • v s = o, x e n si 
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FIG. 3: Abstract setting for the FSI problem considered in the present study. 



p(dtVf + Vf ■ Vvf) =V ■ 07, xettj, 
p(d t v s +v s ■ Wv s ) =V • cr s , xeil s , 

where v denotes the velocity vector, t the time, p the density, and <r the Cauchy stress tensor. 
The no-slip condition is imposed on the fluid-wall boundary, namely 

Vf = V w , x e IV, 



(2) 



(3) 



where Vw denotes the wall velocity. The kinematic and dynamic interactions between the fluid and solid phases are 
determined by continuity in the velocity and in the traction force at the fluid-structure interface, namely, 



v f = v s , xeT T , 



x e Fi, 



(4) 
(5) 



where n denotes the unit normal vector at the interface. 

In the practical numerical procedure based on the full Eulcrian perspective, instead of separately partitioned two 
velocity fields v f and v s respectively in f and in tt s , it is convenient to introduce a monolithic velocity vector v 
applied to the entire domain 0(= r2/(Jf2 s ). In multiphase flow simulations, one set of governing equations for the 
whole flow field, known as a one-fluid formulation [81], is often employed to be discretized on a fixed grid. In the 
present study, such an idea is applied to the fluid-structure system by using v, that is here referred to as a one- 
continuum formulation. The one-continuum formulation would immediately satisfy (4) because v is supposed to be 
continuous across the interface Fj . Following the volume-averaging procedure [72] , we establish the velocity field v as 



V = (1 - (j> s )vf + (j) s V s 

where (f> s is the volume fraction of solid inside a computational cell: 
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(6) 



(7) 



where A denotes the grid size, the suffices x, y and z stand for the respective directions, and I s the indicator function 
defined by 



I s (x) = (l p° r X S 
v ; [0 for x e fi/. 



(8) 



We may regard the volume fraction <p s as a smoothed Heaviside function at the grid scale. The distribution of the 
volume fraction reveals <j) s = in fluid, (f> s = 1 in solid, and < <p s < 1 for the grid involving the fluid-solid interface. 
As explained in §IVB 1, the modulus of its gradient |V0 S | is regarded as a smoothed one-dimensional delta function 
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at the grid scale A x . Namely, for the grid involving only fluid or solid, |V<fc| is zero, whereas for the grid involving 
the interface, |V0 S | reveals a peak of the order of A" 1 . The isoline at <j> s — 1/2 represents the interface Tj (sec figure 
1(b)). In several approaches for multiphase flow simulations (e.g., VOF [23] and level set [6, 70] methods), the viscous 
stress is often written in a mixture form: i.e., the smoothed Heaviside function H smears out the viscosity such as 
/i m i x = (1 — H)/ii + iJ/i2 to remove the discontinuity of the stress across the interface Tt. The mixture representation 
is employed in the present study. For incompressible continua, the pressure p may be regarded as of a Lagrangian 
multiplier imposing the solenoidal condition over the whole velocity field. The Poisson equation will be solved to find 
the pressure field p, written in the one-continuum form, over the entire domain Q. Taking the weighted average with 
respect to <p s , we write the mixture stress er as 

er = ~pl+ (1 - <i>s)o-' f + (j> s cr' s , xe!!, (9) 

where I denotes the unit tensor, and the prime on the second-order tensor stands for the deviatoric tensor, e.g. 
T" = T— |tr(T)J for a tensor T. Since <fi s is smoothed at the grid scale and er is supposed to be smoothly distributed 
over the entire domain, the expression (9) at <f> s = 1/2 would satisfy the continuity of the traction force (5). To advect 
the scalar field <p s . the purely Eulcrian interface method is employed. Throughout the paper, the fluid component is 
assumed to be Newtonian, and thus the deviatoric stress of fluid is given by 

er} = 2n f D', (10) 

where fif denotes the dynamic viscosity of fluid, and D(= |(Vd + Vv T )) the strain rate tensor. Instead of (1) and 
(2) with (4), (5) and (9), we solve the following equations in the one-continuum form over the entire domain: 

V-v = 0, xen, (11) 

p{d t v + v- Vv) = -Vp + V- {2(1-^)^/13' + <j) s cr' s }, x € fl, (12) 

<9 t s + t>-V0 s = 0, x 6 ft. (13) 

The deviatoric stress a' s of solid is dependent on the constitutive law. The incompressible Mooney-Rivlin law [20, 
48, 61] involving a nonlincarity with respect to B (here B = F ■ F T denotes the left Cauchy-Green deformation 
tensor, F = dx/dX the deformation gradient, x the current coordinates, and X the reference coordinates [4]) is 
incorporated into the finite difference method. The constitutive equations and the solution algorithm will be presented 
in the following sections. 

B. Constitutive equations for solid 

In most of finite element computations, the constitutive equations of hyperelastic material are written over the 
reference configuration. The hyperelastic constitutive law is usually provided by using the first or second Piola- 
Kirchhoff stress tensor, which is differentiated with respect to the reference coordinates X in the momentum equation. 
It is suited to the Lagrangian frame. By contrast, in the Eulcrian approach, the equation set is written over the current 
configuration. Therefore, the constitutive equations are needed to be temporally updated on the fixed grid without 
using the Lagrangian mesh. In this section, we describe the constitutive law in the Cauchy stress form, and the 
transport of the left Cauchy-Green deformation tensor field, which are the core elements of the present approach. 

1. Cauchy stress expression of the incompressible Mooney-Rivlin law involving nonlinearity up to 0(B 2 ) 

We consider incompressible visco-hyperelastic materials undergoing only the isochoric motion. The deviatoric 
Cauchy stress of solid is expressed as 

*' s = %l s D' + a' sh , (14) 

where the first term on the right-hand-side corresponds to the viscous contribution with dynamic viscosity (j, s . The 
second term a' sh corresponds to the hyperelastic contribution to be derived below. 
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To formulate the constitutive equation, we refer to the general theories [19, 66, 80] of constrained material. Choosing 
the Mooney-Rivlin expression [48, 61], and considering the nonlinearity up to 0(B 2 ) in the deviatoric Cauchy stress, 
we write the hyperelastic strain energy potential W as 

W = ci(T c - 3) + c 2 (5 c - 3) + c 3 (T c - 3) 2 , (15) 

where lc and He denote reduced invariants [20] of the right Cauchy-Green deformation tensor C — F T ■ F defined 
by 

l C = — TJo, n C = — TH, (16) 



where 



I c =tr(C), II c = ^{I^-tr(C-C)}, III C = det(C). 



Utilizing the equivalence of the invariants between the left and right Cauchy-Green deformation tensors [24], we write 
the deviatoric Cauchy stress tensor as 

~ - -F {2- — --^ + 2- ^-£- + 2- ——3— \ • F 



sh det(F)" \~dl c dC ' ~dll c dC ' ~ dlll c dC J * (17) 
= (2 Cl B + 2c 2 (tr(B)B B B) + 4c 3 (tr(B) - 3)B)'. 

We will give several demonstrations afterward for some specific cases based on the linear Mooney-Rivlin, neo-Hookean, 
and incompressible Saint Venant-Kirchhoff materials. Note that all these materials obey (17). For linear Mooney- 
Rivlin material (03 = 0) [24, 48, 61], which is often used as the biological material, Cauchy stress becomes 

a' sh = 2ciB' + 2c 2 (tr(B)B - B ■ B)' . (18) 

The neo-Hookean material is a particular case of the linear Mooney-Rivlin material with the coefficients 

G 

ci = — , c 2 = c 3 = 0, 

where G denotes the modulus of transverse elasticity. The corresponding Cauchy stress is 



= GB'. (19) 



As another typical hyperelastic material, we consider Saint Venant-Kirchhoff material [4, 13, 63], which often models 
a thin but finite volume membrane. The constitutive equation is expressed as a simple extension of Hooke's law, as 
defined by 

s = K^ME)i + ii<~.. .., k. 

where S is the second Piola-Kirchhoff stress, X^^, are the Lame constants, and E = |(C — J) is the Green- 

Lagrange strain tensor. Although the Saint Venant-Kirchhoff material is usually referred to as dilatable, we implement 
the incompressible Saint Venant-Kirchhoff model, as defined in [18]. The deviatoric Cauchy stress of the incompressible 
Saint Venant-Kirchhoff material is expressed as 

' = %i( tr (B) - 3)fl' + ^(B ■ B B)'. (20) 



det(F) 



Substituting the relations 



2 ' 8 

into (17), we can readily realize that the constitutive law (20) falls within the class of nonlinear Mooney-Rivlin laws. 
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2. Transport of left Cauchy- Green deformation tensor field in the Eulerian frame 



Once the coefficients c±, ci and C3 are given, the constitutive equation (17) is expressed as a function of the left 
Cauchy-Green deformation tensor B. If the tensor field B is determined in a purely Eulerian manner, all the equations 
will be closed in the Eulerian form. Utilizing the fact that the upper convected time derivative of B is identically 
zero [4], one may use the following transport equation to update the B field: 

d t B + v VJ3 = L B + B L T , (21) 

where L(= (Vv) T ) denotes the velocity gradient tensor. For the initially unstressed solid, the initial condition is 
given by B — I. It should be noticed that, however, it is quite cumbersome to solve (21) from a numerical viewpoint, 
because B exhibits rough distribution in the fluid domain Qf [46]. The fluid element subject to a shearing motion 
is likely to elongate toward the extensional direction. Such an elongation causes a temporally exponential growth 
of some components of B. To avoid the numerical instability brought by the exponential growth, instead of B, we 
update the modified left Cauchy-Green deformation tensor B = <fi^B (a > 0), which yields 

d t B + v VJ5 = L B + B L T ', (22) 

with the initial condition B = 4>"I. Because of B = for (j> s = 0, the numerical instability is evaded in the fluid 
domain ilf. In the hyperelastic constitutive law (18), we consider the term up to 0(B 2 ), of which the contribution to 
the phase averaged stress 4> s (Ts in (9) is proportional to <fi s B 2 = <\>\r 2a B 2 . In order to avoid division by zero in the 
fluid region <p s = 0, the exponent 1 — 2a on <f> s should be non-negative, indicating a < 1/2. In the present study, we 
choose the largest value a = 1/2. Further, to avoid the inevitable exponential growth at the cell near the interface F/ 
containing the fluid-solid mixture, and to obtain a viable compromise between the numerical consistency and stability, 
we introduce a threshold </> m i n and enforce B = where <f> s < m i n . In the present study, we set (j) m i n between 0.01 
and 0.1. 



3. Characteristics of the present approach in treating the solid stress 

From (17), the resulting deviatoric stress of solid multiplied by <f> s , which is involved in (12), is expressed as 

</>.</„ =2^ sf x s D' + (2 Cl <\>\l 2 B 

+ 2c 2 (tr(B)B ~B-B) + 4c 3 (tr(B) - 30 5 1/2 )B)', 

which can be evaluated together with the temporally updated B from (22). 

It should be noted that when the material points are tracked in the Lagrangian way, the relation (21) is identically 
satisfied via the change in the relative position of the adjacent material points (see figure 2(a)). Thus, in the pure 
Lagrangian approach, it is not necessary to introduce a transport equation describing the deformation level such as 
B. On the other hand, in the Eulerian approach, which does not rely on the material point, we must introduce 
the deformation level on the fixed mesh. The present approach is characterized by the introduction of the transport 
equation of B, which is temporally updated in the Eulerian frame (see figure 2(b)). 

III. SIMULATION METHODS 
A. Overview 

Once the initial field of the solid volume fraction (j> S Q is given over the entire domain f2, the present Eulerian method 

enables one to carry out the FSI simulation without mesh generation procedure. In order to prescribe the (f> S Q field, 

it is only required to numerically compute the ratio of the occupied solid to each control volume from the initial 

geometry as a preproccss. If the system is initially at rest and unstressed, one can launch the simulation by setting 

the initial conditions of the velocity vector, pressure, and modified left Cauchy-Green deformation tensor to v = 0, 
1/2 

p = 0, and B = q>s I, respectively. 

Instead of the deviatoric stress er' (with a prime), we may use the following pseudo stress to make some discretized 
expressions for the individual stress components simple: 

& = (p f + {/is - n f )<f, a ) (V« + Vu T ) + cf> s fr sh , (24) 
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where 



4>„a Bh = 2c 1( /)l /2 B + 2c 2 (tr(B)B - B ■ B) + 4c 3 (tr(B) - 3</> s 1/2 )B. 
Introducing a pseudo pressure p, instead of (12), we solve 

p (d t v + v ■ Vu) = —Vp + V • cr. 



(25) 



(26) 



The actual pressure and deviatoric stress read 



p = p 



tr(cr) 



cr = cr 



tr(cr) 



Hereafter, the formulation will be given on the two-dimensional (2D) system. The extension to the three-dimensional 
(3D) system is straightforward [69]. The basic equations (11), (12), (13), (22) and (23) are solved by a finite difference 
method on a fixed Cartesian grid. The quantities v, p, <p s , and B are temporally updated in the Eulcrian frame. 
The entire domain is discretized by a uniform square mesh. We follow a conventional staggered Marker- And-Cell 
(MAC) cell arrangement [21], where the velocity component is located on the cell face, and the pressure at the cell 
center (see figure 4(a)). The symmetry in the tensor B and the two-dimensionality of the system imply B zz = 



B xz — B zx 

and B xy . The diagonal components of B and the solid volume fraction 
non-diagonal components of B are on the cell apex (see figure 4(b)). 



B yz = B zy = 0, and B yx = B xy . Hence, (22) is solved for the three independent components B xx , B yy 

are defined on the cell center, while the 





FIG. 4: Schematic figure of computational grid with the mesh size of A x x A y , here A^ = A y . (a) left panel: Definition points 
of the velocity components v x , v y , and the pressure p. (b) right panel: Definition points of the solid volume fraction <f> B , the 
stress components a xx , a yy , o xy , and the modified left Cauchy-Green deformation components B xx , B yy , B xy . 

As opposed to usual methods of computational structure dynamics, for instance, the GMRES approach [62] and 
the weak compressibility approach [7, 28], the pressure Poisson equation is solved to exactly satisfy the solcnoidal 
condition (11) over the entire domain f2 likewise widely-used incompressible fluid flow algorithms. In solving the 
discretized Poisson equation, the fast Fourier transform is used to ensure high accuracy as well as high efficiency. 



B. Time-stepping algorithm 

Here, the time-stepping algorithm to update the variables at the (n + l)-th time level from the n-th time level 
is briefly explained. Following the Simplified MAC method [1], corresponding to a standard incompressible fluid 
flow algorithm, with an incremental pressure correction applied to the finite difference scheme, we decompose the 
time-stepping into three steps. 

In the first step, the second-order Adams-Bashforth scheme [5] is applied to explicitly updating the volume fraction 
cj>i n+1 ^ and the modified left Cauchy-Green deformation tensor B( n+1 ); 

(n+l) = (n) _ (At) Q„(n) . ^(n) _ ^(n-l) . , (27) 
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£?(»+!) = £?(») - (At)||«W ■ V£? (n) - \v {n ~ 1] ' VB (n_1) 



| (-£(") • f") - flW • L T (")) (28) 

1 

2 



where (At) denotes the time increment, and the superscript (n) stands for the n-th time level. It should be noted 
that (At) is chosen such that the Courant-Friedrichs-Lewy (CFL) condition is satisfied i.e. the CFL number based 
on the larger velocity scale between the maximum advection speed maxdw^j, \v y \) and the linear elastic wave speed 
\j2(c\ + C2)/ p is 0.1 or less for all the present computations. 

In the second step, the second-order Adams-Bashforth and Crank-Nicolson schemes [5] are applied to iteratively 
calculating the unprojected velocity vector v* and stress tensor <x*: 



v * = v (n) _ 



(At)Vp(™) 



- (At) Q»W • W"> - \v {n - 1] • W"" 1 ^ 



(29) 



a* = (p f + (p s -p s )^ +1) ) {Vv*+Vv* T ) + (<f> s <T sh y n+1 \ (30) 

where the solid stress is given by (25). 

Finally, pressure, solenoidal velocity vector, and stress tensor are updated during the projection step: 

p(n+l) = £(«) + ^ ( 31 ) 

„(n+i) = v * _ (At)V(p, (32) 

= & * _ 2(At)(p f + (p s - M/ )0(" +1 ))(VV^), (33) 
where the incremental pressure if is determined by solving the Poisson equation 

W=^f (34, 



C. Spatial discretizations 



The spatial derivatives arc approximated by the second-order central differences, except for those of the advection 
terms in (13) and (22), to which the fifth-order WENO scheme [36, 45] is applied. For the momentum equation, 
following the spirit in [34, 38], we discretize the advection terms to satisfy the identity V • (vv) = (v ■ V)u + v(V ■ v) 
in the discretized space, that would make the energy highly conserved. The finite difference descriptions are detailed 
in Appendix A. 

We here focus on the discretization of the right-hand-side of (22), that is important to accurately describe the 
isochoric solid deformation. We take care of the difference between the definition points of the diagonal B xx , B yy and 
non-diagonal B xy components as illustrated in figure 4(b). The incompressibility det(JB) = 1 implies B xx B yy — B 2 y = 
(j) 2a . Since we choose a = 1/2, we find 

A J£d 2 x (B xx B yy - B 2 xy ) = A JJ^ X ^ = 0, (35) 

In consideration of the time derivative of each component of B in (22), the left-hand-side of (35) is given by 

d 2 x ^(d t B xx )B yy + B xx (d t B yy ) — 2(d t B xy )B xy ^ 

d 2 x {(-v ■ VB xx )B yy + B xx (-v ■ VB yy ) - 2(-v ■ VB xy )B xy } 
a (36) 



d 2 x {e x (L B + B ■ L T ) ■ e x B yy + e y ■ (L ■ B + B ■ L T ) ■ e y B xx 
2e x ■ (L-B + B-L T )-e y B xy }, 
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where e denotes the unit vector. Applying Gauss' divergence theorem together with the kinematic condition n ■ v = 
at the rigid wall T w and the solenoidal condition (11), we rewrite the first term in the right-hand-side of (36) as 



I dec n-v{B xx B yy - B 2 } = 0. 



Hence, in view of the solid volume conservation, the following relation should be satisfied: 

d 2 x {e x (L B + B L T ) ■ e x B yy + e y ■ (L ■ B + B ■ L T ) ■ e y B xx 
n (37) 

- 2e x (L B + B L T ) ■ e y B xy ) = 0. 

We choose the interpolation and the finite difference formulae to satisfy the integral relation (37) in a discrete form. 
For a quantity q, let us introduce finite difference operators 6i and Sj, of which the indices i and j correspond to 
discrctized coordinates along the respective directions x and y, such as 

Si(q)\i,j = q i+ i d - Ui-ij, Sj(q)U,j = - ( 38 ) 

and a four-point interpolation operator denoted by double ovcrline such as 

The integral J„ d 2 x fg provides the equality in a finite volume representation 



2 2 ^ 2 2 ' 

i j i j 



where the quantities / and g are defined at the cell center and at the cell apex, respectively, and assumed to vanish 
at the boundary Tw We write each component of L ■ B + B ■ L T involved in (37) as 



where 



(e x (L B + B L T ) ■ e x ^j — 2L xx jjB xx jj + 2L xy B xy | . . , (41) 
(e y ■ (L ■ B + B ■ L T ) ■ e y ^j_ _ = 2L yy ^B yy ^ + 2L yx B xy \.., (42) 
(e x ■ (L ■ B + B ■ L T ) ■ e y ). +hj+ = +L^\ i+hj+ i) B xy ,i+U+\ 

+ L xy,i+^j+^ B yy\ l+ i J+ i + L yx,i+^,j+^ B xx\ l+ l J+ l, (43) 
Si(v x )\i,j T _ &j{v y )\ij 



_ Sj(v x )\ l+ i d+ i S l (v y )\ i+ i J+ i 

J-/^.,, ,il o^I ". , „•_!_! 



*y,i+2>3+2 A ' ^.i+2J+2 A 

Substituting (41)-(43) into (37) together with the relation (40), and considering the solenoidal condition (11), we 
confirm that the requirement (37) would be fulfilled in a discrete form as 

2 y ] y ] A x A y {L xx .ij + L yy ,ij) ( B xx ^jB yy ^j — B^ylij J = 0. 

J =0 
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IV. VALIDATION TESTS 



For the purpose of addressing several computational issues involved in the procedure advocated in §11 and §111, 
firstly, we consider a one-dimensional problem for the motion of the oscillatory parallel fluid-solid layers. Though this 
system is a simple example, it includes the fundamental aspects of the FSI problem and the system allows estimate of 
the numerical accuracy in the present full Eulerian approach as it has analytical solutions. In the second example for 
validation, we make comparisons with the available data of the deformablc solid motion in a lid-driven cavity [97] and 
the particle-particle interaction in a Couctte flow [12]. Also, in the third example, a response of hypcrelastic material 
to the external shear strain is examined to check the reversibility in shape as this aspect is vital for a full Eulerian 
formulation. 



v w (0 = Vw sin mt = Im ( Vw exp(itftf) j , 


Fluid 




► x solid 


2L S 


Fluid 


Lf 



-v w (t) 



FIG. 5: Schematic figure of the parallel fluid-structure layers. The upper and lower plates move in opposite directions with 
temporally sinusoidal velocities to drive the fluid and solid motions. 



A. Oscillatory response in parallel layers of fluid and solid 



As schematically illustrated in figure 5, we deal with the interactive motions of the three (fluid-solid- fluid) parallel 
layers bounded with two oscillatory plane walls. Here, making comparisons with accurate solutions obtained by 
means of a sharp interface approach, we examine the validation and verification of the present Eulerian approach, 
accompanied with a diffuse interface. Supposing homogeneity in x direction, we may omit the x-dependencc of any 
quantity in the theoretical analysis. In the numerical simulation, the periodic condition is applied in x direction. We 
here treat pure hyperelastic material, i.e., (i s = 0. The relations between the velocity v, the displacement u and the 
shear stress a are given by 

d t v = d y a, (44) 

d t u = v, (45) 

H f dyV for fluid (L s < \y\ < L s + Lf), . . 

2(ci + c 2 )d y u + Ac 3 (d y u) 3 for solid (0 < \y\ < L s ), yw > 

with no-slip condition at the upper and lower plates (y = ±(L/ + L s )) 

V w sinwt at y = L f + L s , . , 

— Viysinwf at y = — (Lf + L s ). 
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FIG. 6: The velocity profile in the upper half region of the parallel fluid-structure layers under the conditions of /i/ = 1, 
L s = Lf = 0.5, u) = 7r, and V w = 1. Left panels: the neo-Hookean material with the modulus of transverse elasticity of G = 5. 
Right panels: the incompressible Saint Venant-Kirchhoff material with the Lame constants of AL amd = 7.5 and ^, m j = 5. The 
upper and lower panels show the results at the temporal phases of t%2 = and t%2 — 1.8, respectively. The dashed and solid 
curves respectively represent the linear and nonlinear solutions with the sharp interface, which are determined by means of the 
spectral approach (see Appendix B). The symbols correspond to the present Eulerian simulation results for various number of 
grid points (N x X N y = 8 X 32, 8 x 128, 8 x 512). 



The solid stress expression (46) indicates that the system involving a linear Mooney-Rivlin material with C3 = is 
linear with respect to the displacement u, while that involving the Saint Venant-Kirchhoff material with C3 =/= is 
nonlinear. 

The simulation results based on the present full Eulerian approach are compared with the analytical solution 
obtained by means of the sharp interface approach (see Appendix B for detail). We fix the conditions p = 1, jif = 1, 
L x = 8, L s = Lf = 0.5, lu = 7T, and V w = 1, and vary the number of grid points (N x x N y = 8 x 32, 8 x 128, 8 x 512). 
Initially, the system is at rest. The total computational period is t = 40, corresponding to 20 cycles, and the sampling 
is performed within the last one cycle i.e., t £ (38,40]. Figure 6 shows the velocity profiles for the neo-Hookean 
model (G = 5, i.e., c\ = 2.5, 02 = 03 = 0) and for the Saint Venant-Kirchhoff one (AL amd = 7.5, /i£ amiS = 5, i.e., 
ci = 5.0, C2 = —2.5, C3 = 2.1875). The present simulation results converge to the sharp interface solutions with the 
higher spatial resolution. Figure 6(d) shows the obvious difference in profile between the linear (03 = 0) and nonlinear 
( c 3 7^ 0) solutions at the phase of t%2 =1.8 (here % stands for the remainder). The simulation results in figure 6(d) 
clearly get closer to the nonlinear solution with increasing the number of grid points, indicating that the nonlinearity 
in the solid constitutive law is reasonably captured in the present Eulerian approach. 

The accuracy in the fluid-structure coupling is quantified by the errors in L2 and Loo norms, which are respectively 
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FIG. 7: Upper panels: The error of the velocity in L2 norm versus the number N y of grid points in the vertical (y) direction 
for various temporal phases in the parallel fluid-structure layers. The error is determined from the difference between the 
results of the present Eulerian and sharp interface methods based on (48). Lower panels: The absolute slope in the plot of 
the error versus N y . The local slope is determined from (49). The left and right panels correspond to the neo-Hookean and 
incompressible Saint- Venant Kirchhoff materials, respectively. The conditions are the same as those of figure 6. 

denned as 

n^ii(iv 2/ ) = |^E(^ ) -^ ) ) 2 

\ N vU 

where Vj and , respectively, denote the present result and the sharp interface solution on the node yj = {j—\)/N v . 
The log- log plots of the L2 and Loo errors versus the number of grid points are shown in figure 7(a)(b) and figure 
8(a) (b), respectively. The local slopes are obtained therefrom, and shown in figure 7(c) (d) and figure 8(c) (d). Here, 
the local slopes are determined from the following approximations 

aiog(|[L 2 ||) log (\\L 2 \\(V2N y )) - log {\\L 2 \\{N y /V2)) 

d\og{N v ) 1 y) ~ log2 ' 

aiogfllLooll) ^ iog(||£ 00 ||(v^iv,))-iog(||i:oo||(iv,/V2)) 

d\og{N y ) 1 v >~ log 2 

The slope indicates the degree of the accuracy in the fluid-structure coupling. In both the L 2 error (corresponding to 
a global indicator) and the L^ error (corresponding to a local maximum indicator) are nearly proportional to -/V" 1 . 
Since the second-order finite difference is applied to describing the spatial derivatives, this near-first order trend must 
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FIG. 8: Same as figure 7, but in norm. 



have resulted from the mixture stress expression (9) involving the first-order accuracy locally at the interface, which 
dominates the global degree of accuracy. Note that we also investigated the grid convergence of the shear stress, and 
confirmed that the accuracy at the interface is of the first-order with respect to the grid size. 

Figure 9 shows the sensitivity of the wall friction amplitude T WtIlas on moduli. We fix the conditions ptf = 1, /i s = 0, 
L a = Lf = 0.5, lj = 7r, and V w = 1. Likewise the computations in figure 6, the total computational period is set to 
20 cycles. The root-mean-square of t w was sampled over the last one cycle. The results of the linear Mooney-Rivlin 
model with (02 = ci/2, c-s = 0) as well as the neo-Hookean (c2 = C3 = 0) and Saint Venant-Kirchhoff (02 = — ci/2, 
C3 = 7ci/16) models are plotted as a function of 2(ci + cq). Because the deformed motion of solid behaves like as 
the spring-mass system, the plot of the wall friction amplitude versus 2(ci +C2) reveals the non- monotonous resonant 
behavior. As long as the solid strain is sufficiently small, the nonlinearity involved in the constitutive law is negligible, 
and therefore the linear assumption is justified. Indeed, the curve of the nonlinear solution approaches the linear 
solution with increasing 2(ci + C2) since the solid strain is suppressed for the stiffer material. By contrast, for the 
smaller 2(ci +C2), the discrepancy between the linear and nonlinear solutions becomes more obvious. It is because the 
larger strain makes the nonlinear system effectively stiffer as implied by (46). All the results of the present Eulerian 
approach are in good agreement with the sharp interface solution. The present approach is confirmed to capture the 
resonance behavior resulting from the dynamic interaction between the fluid and solid motions, and the nonlinearity 
in the solid constitutive law. 



B. Comparison with independently conducted FSI analyses 

We here make comparisons with two well-validated FSI analyses. In the constitutive law for (visco-)hyperelastic 
material, one has [i s = [if, c\ 7^ and ci = C3 = 0, and the other has fi s = 0, C2 ^ and c\ = C3 = 0. 
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FIG. 9: Response curve of the skin friction amplitude T w ,rms m the motion of the parallel fluid-structure layers as a function of 
2(ci +C2) with [if = 1, L s = Lf = 0.5, ui — n, and V w = 1. The symbols correspond to the results based on the present Eulerian 
method with the 8 x 256 mesh. The circles, crosses, and squares represent results of the neo-Hookean (02 = C3 = 0), linear 
Mooney-Rivlin {c% = ci/2, C3 = 0) and incompressible Saint Venant-Kirchhoff (02 = — ci/2, C3 = 7ci/16) models, respectively. 
The curves represent the sharp interface solutions, which are determined by means of the spectral methods (see Appendix B). 
The dashed curve represents the linear solution for the neo-Hookean and linear Mooney-Rivlin materials with C3 = 0, while the 
solid curve the nonlinear solution for the above-mentioned Saint Venant-Kirchhoff material. 



1. A solid motion in a lid-driven cavity flow 



We perform full Eulerian simulations of deformable solid motion in a lid-driven cavity with the same setup and con- 
ditions as Zhao et al. [97], who employed mixed Lagrangian and Eulerian approach. The initial setup is schematically 
illustrated in figure 10(a). The size of the cavity is L x x L y = 1 x 1. Initially, the system is at rest. The unstressed 
solid shape is circular with a radius of 0.2, and centered at (0.6, 0.5). At t = 0, to drive the fluid and solid motions, the 
top wall starts to move at a speed of Vw = 1 in x direction. The no-slip condition is imposed on the walls. The solid 
component is neo-Hookean material. The material properties are p = 1, fj,f = jj, s = 10~ 2 , c\ = 0.05 and C2 = C3 = 0. 

Figure 10 visualizes the particle deformation and the flow field for eight consecutive time instants. The dashed 
curve in figure 10 represents the outline of the particle obtained by Zhao et al. [97], in which they computed the 
solid deformation on the Lagrangian mesh. The solid lines represent the instantaneous particle shapes, corresponding 
to the isoline at <fi s = 1/2, obtained by the present full Eulerian simulation. The dotted material points are tracked 
just to transfer images of the particle deformation, but we did not use these material points for computing solid 
stress and strain. The particle moves and deforms driven by the fluid flow, and exhibits highly deformed shape when 
the particle approaches the wall. It should be noticed that no special artifact for avoiding a particle-wall overlap is 
implemented into the present method because the particle- wall hydrodynamic repulsion is likely to be brought by the 
soft lubrication effect [67] due to the geometry change via the particle deformation. The solid shapes obtained by the 
present Eulerian simulation are in excellent agreement with the well- validated result [97]. 

In addition to the comparative study, we address grid convergence issues below. We trace the centroid x c = (x c , y v ) 
of the particle, which is evaluated from the approximation 

N X Ny 

y^y^ AsA^ x it j <f> s (x itj ,t) 

x c (t) w 1=1 j= \ . (50) 

i=l 3=1 

Figure 11 shows the trajectory of the centroid x c in a time range of t £ [0,20] for various number of grid points 
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FIG. 10: Comparison of the solid deformation in the lid-driven flow with the simulation result [97]. The dashed outline represents 
the result of Zhao et al. [97], in which the Lagrangian tracking approach was employed to describe the solid deformation. The 
solid outline, the dotted material points and the streamlines correspond to the present simulation results based on the full 
Eulerian approach with a mesh 1024 x 1024. 



(N x x N y = 64 x 64, 128 x 128, 256 x 256, 512 x 512, 1024 x 1024). The trajectories clearly exhibit a convergent trend 
to the curve of the highest spatial resolution. 

To quantify the grid convergence behavior, the distance of the particle centroid x c with respect to that of the 
highest resolution N x = 1024 is monitored. We evaluate the errors in L2 and norms respectively from 



\\L 2 \\(N X ) =|i jTdi \x c (t,N x )-x c (t,N x = W24)f 

\\L 00 \\(N X )= max |as c (i, N x ) - x c (t, N x = 1024)|. 
te[o,T] 



(51) 



Figure 12 shows the L 2 and errors as a function of N x . Both the errors are nearly proportional to N~ x , indicating 
the first-order accuracy. The first-order accuracy involved in the fluid-structure coupling as described in §IVA is 
reflected on the slopes in these plots. 

In addition to the particle centroid, to further examine the grid convergence behavior, we here perform modal 
analyses of the particle deformation. The distance from the particle centroid x c to the interface xj = (xi,yi) is 
written as 

R(6) = \x T - x c \, (52) 
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FIG. 11: Trajectories of the solid centroid in the lid-driven flow in a time range t £ [0, 20] for various number of grid points. 
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FIG. 12: The errors of the particle centroid in L2 norm and in Loo norm versus the number N x of grid points in the lid-driven 
flow. 



where 9 is found to satisfy the relations 



xi -x c yi - y c 

■ r , Sin = : 



\Xi - x c 



\Xi - x c 



The distance is written in a Fourier scries form 



R{9) = Rq + (R cn cos n9 + R sn sin n9) , 



(53) 



(54) 



where R n denotes the n-th order deformation mode. The deformation modes arc uniquely determined via the orthog- 
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onality in the cosine and sine functions from definite integrals 

i r 2 " 



R cn =- d6 R(6) cos nG, R sn = - dO R(9) sin nO. 

7T Jo K JO 




FIG. 13: Time history of the n-th order modal amplitude \Rn\ of the particle deformation in the lid-driven cavity flow for 
various number of grid points, (a) n = 0, (b) n = 1, (c) n = 2, and (d) n — 3. 

In the present Eulcrian approach, the fluid and solid phases are distinguished by the solid volume fraction (p s , and 
there is no explicit quantity describing the angular profile of R. Instead of the circumferential integral in (55), we 
will apply the area integral to evaluating the deformation mode. Let us consider the following relation for a function 

2 \eR{e)f{9)= f[d 2 x6(\x- Xl \) f(6), (56) 

J JQ 

where <5 stands for a one-dimensional Dirac's delta function, which is related to the gradient of the solid indicator 
function / s , namely, 

VI S = — nS(\x — xi\), 
where n denotes the unit normal vector pointing towards the fluid and is given by 

VI S 

n WIS 
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Hence, the delta function used in (56) is expressed as 



<y(|sc-sej|) = |VJ a |. 

Implementing (56) in the finite difference approach, one must smooth the delta function at the grid scale. Upon using 
the solid volume fraction (f> s , which is regarded as the solid indicator function smoothed at the grid scale, we obtain 
the approximation 

tfflas-as/DwIV^I, (57) 
and then evaluate the definite integrals in the summation form 



f cW R(9)f(6)^f^A x A y \V ( l> s , i , j \f(e i , j ), 
J ° i=l ,=i 



(58) 



together with (50) to find Oij = 9(xj,ij, x c ). For n > 1, we write the modal amplitude as |i? n | = -J R? cn + R? sn . 

Figure 13 shows the temporal evolutions of the modal amplitude \R n \ (n = 0, 1, 2 and 3) of the particle deformation 
for various number of grid points. The largest elongation (n — 2 mode) of the particle is observed at about t = 5 
when the particle is in the proximity of the moving wall, and the synchronized increases in \R n \ are found for different 
modes. With increasing the number of grid points, the profiles for each mode n settle to corresponding convergent 
curves, suggesting the verification with respect to the deformation of the particles. 




FIG. 14: The errors of the n-th order modal amplitude \R n \ (n = 0, 1, 2, 3, 4 and 5) of the particle deformation (a) in L2 norm 
and (b) in norm versus the number N x of grid points in the lid-driven flow. 



In a similar manner to (51), the L2 and Loo errors of the modal amplitude with respect to that of the highest 
resolution N x — 1024 arc quantified as follows: 

\\L 2 \\(N X ) = I At \\Rn\(t,N x ) - \Rn\(t,N x = 1024) | 2 1 , 

\TJo J (59) 

\\Loo\\(N x ) = max \\R n \(t, N x ) - {R^t, N x = 1024)1. 

te[o,T]' 1 

Figure 14 shows the L2 and Lqo errors as a function of N x . Again, both the errors are nearly proportional to iV" 1 , 
indicating the first-order accuracy for capturing the particle deformation. 



2. Two particles interaction in a Couette flow 



We here make a comparison with the available numerical analysis of the interaction between two dcformablc particles 
in a Couette flow performed by Gao & Hu [12], who adopted body-fit Lagrangian mesh. The computational extent 
is L x x L y = 8 x 4, which is the same as [12]. Initially, the system is at rest. Two unstressed solid particles are 
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FIG. 15: Comparison of particle-particle interactions in a Couette flow with the simulation result [12]. The dashed outline 
represents the result of Gao & Hu [12], in which the body-fit Lagrangian mesh was used to solve the FSI problem. The dotted 
material points and the solid outline correspond to the present simulation results based on the full Eulerian approach with a 
mesh 1024 x 512. At t — 0.8, initial deformation; At t = 4.0 and t = 5.6, "roll over" interacting mode; At t = 33.2, "bounce 
back" interacting mode [12]. 



initially circular with a radius of 0.5, and centered at x c ^a = (2,2.5) and x c .b = (6, 1.5) as depicted in figure 15(a). 
The upper and lower plates located at y = 4 and y = 0, respectively, start to move impulsively to drive the fluid and 
solid motions at speeds of V^ ppcl = 1 and V^ wor = — 1 in x direction. 

The no-slip condition is imposed on the plates, while the periodic condition is applied in x direction. The solid 
component is purely hyperelastic. The material properties arc p = 1, /x/ = 20, fj, a = 0, C2 = 40 and C\ — c 3 = 0. 

Figure 15 visualizes the two-particle shape for five time instants. The dotted markers are, again, to represent the 
solid deformations and those markers are not used for computing solid stress or strain. The arrows at the particle 
centers are the instantaneous translating velocity vectors. The dashed curve in figure 15 represents the outline of the 
particles obtained by Gao & Hu [12]. The particles experience somehow complicated interactions involving the "roll 
over" and "bounce back" modes as examined in [12]. The solid shape obtained by the present Eulerian simulation is 
again in agreement with the well- validated result [12], indicating that the particle-particle interaction is also reasonably 
captured by the present approach. 

Figure 16 shows the temporal evolution of the y c -position of the particle centroid, which is evaluated from (50), for 
various grid resolutions (N x x N y = 128 x 64, 256 x 128, 512 x 256, 1024 x 512). In the full Lagrangian computation 
[12], the finite element mesh is refined within the particle-particle gap, whereas in the present Eulerian simulation, 
the grid size is uniform and fixed. When the plot shows peaks around t = 3.0, t = 16.0 and t = 20.0, the gap between 
the particles is narrow, and the particle undergoes relatively strong hydrodynamic force owing to a squeezing effect. 
Such a narrow-gap effect is less resolved by the present method than the full Lagrangian method especially for the 
low spatial resolution cases, that is reflected on the larger deviations from the result by Gao & Hu [12] preferentially 
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FIG. 16: Variations of particle imposition as functions of time for various number of grid points. Comparison with the result 
of Gao & Hu[12]. 



at the peaks. In the higher spatial resolution, the profiles of the present simulation get closer to the full Lagrangian 
result [12]. 

C. Reversibility in shape of hyperelastic material 

The hyperelastic material generally exhibits reversibility in shape when it is released from stress. In the total 
Lagrangian method using the finite element mesh, since the tracked material point links both the reference and 
current configurations, the reversibility can be captured with little difficulty. By contrast, the Eulerian fixed grid 
point retains no information on the reference configuration. Therefore, one may raise a shortcoming that the Eulerian 
approach is likely to lose the information about the original shape once the material is stressed to deform. We here 
perform a reversibility examination. 

1. For a circular particle 

We here deal with a shear flow between two plane plates involving a hyperelastic particle. The distance between 
the plates is L y = 2. The computational extent in x direction is set to L x = 8. The upper and lower plates are 
located at y = 1 and y = — 1, respectively. Initially, the system is at rest. An unstressed solid particle is initially 
circular with a radius of 0.75, and centered at the middle position (0, 0) between the plates as depicted in figure 17(a). 
The no-slip condition is imposed on the plates, whereas the periodic condition is applied in x direction. We fix the 
material properties p = 1, [if = 1 and (i s = 0. We consider two kinds of materials: one is the linear Mooney-Rivlin 
material with c\ = 4, ci = 2 and C3 = 0, and the other is the Saint Venant-Kirchhoff material with A£ am(l = 6 and 
f*Lm* = 4 (i-e-, ci = 4, c 2 = -2 and c 3 = 1.75). 

The system motion is controlled as follows. Within a period of < t < 4, the upper and lower plates move at 
speeds of V^ ppcr = 1 and V^ wcr = — 1 in x direction, respectively, to drive the fluid and solid motions. After t = 4, 
the moving plates stop (i.e. Vj^ ppor = V^ wer = 0) to release the particle from the shearing force. 

Figure 17 visualizes the particle deformation and the flow field for six consecutive time instants. As the shear flow 
is induced by the moving plates, the shearing force is imposed on the solid particle, and causes the particle elongation 
toward the extensional direction. In the transient state during the development of the deformation, it is observed 
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FIG. 17: Snapshots of the velocity (arrows) and vorticity (color) fields involving a circular particle in the imposing-releasing 
shear flow between two parallel plates. The number of grid points is 1024 x 256. The upper and lower plates move at speed of 
1 and —1 within a period of t G [0, 4], and then stop after t = 4. The solid obeys an incompressible Saint Venant-Kirchhoff law 
with ix a = 0, ALmd = 6 and /il ala6 = 4. 



in figure 17(b) (c) that the transverse elastic waves travel inside the solid, and are reflected by the fluid-structure 
interface. The wave amplitude is damped through the repetitious reflections with time as shown in figure 17(d)(e). 
As examined in [12], the elastic wave propagation inside the particle may play an important role on the deformation. 
As shown in figure 17(e), the vorticity inside the particle at t = 4 is negative, indicating that the particle experiences 
a tank-treading like motion. After the shearing force is released by setting the wall velocities to be zero at t = 4, 
the fluid flow rapidly decays and the deformed particle gradually recovers the circular shape. At t = 6 as shown in 
figure 17(f), the vorticity in the bulk fluid is almost zero, while the non-zero vorticity forms near the fluid-structure 
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interface, indicating the particle shape is under recovery. 




FIG. 18: The budget of the kinetic-energy transport (60) in the imposing-releasing shear flow. The conditions are the same 
as those of figure 17. The dotted, solid, dashed, and dashed-dotted curves correspond to the energy input rate I, the strain 
energy rate — e s , the energy dissipation rate — £/, and the kinetic-energy transport — dE/dt, respectively. The each component 
is provided in (61). The dashed-double-dotted curve corresponds to the summation of the left-hand-side terms of (60). 

The complete recovery in the particle shape is established when the byperelastic strain energy potential returns to 
the initial state. Therefore, the energy transfer between the fluid and solid phases is relevant to it. We here check 
whether the energy transport is numerically conserved during the simulation. In the system addressed in this section, 
the budget of the kinetic-energy transport is written as 

AE , s 

Z-e s -e f - — =0, (60) 

where X, s s , ef, and E, respectively, denote averaged quantities of the energy input rate, the strain energy rate, the 
energy dissipation rate, and the kinetic-energy, expressed as 



Ly \ " W \ dy / r ..ppcr w \ dy / 

e s = (4> a D' : On . (61) 
e f =2/i / ((l-^)r>':£)') a , 

E =f <« ■ o) n , 

where (■■■)r w stands for the average over the wall, and (...)f2 for the average over the entire domain ft. Figure 18 
shows the time history of the each contribution in the left-hand-side of (60). As the flow evolves, the particle deforms, 
and thereby the particle stores the strain energy as indicated by the solid curve with negative value in figure 18. 
After the walls stop at t = 4, the particle releases the strain energy. The double-chained curve in figure 18 shows 
the summation of the left-hand-side terms of (60). Its absolute value, corresponding to the numerical error, is less 
than 10 , which is much smaller than the variation of the contributions of the individual terms. The system is well 
conserved during the simulation in view of the energy balance, because the equality of (60) is almost fulfilled. It is 
important to emphasize that the numerical energy conservation hinges upon the finite difference schemes, and the 
method proposed by Kajishima [38] is employed in the present study. 
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FIG. 19: Material point distribution in the imposing-releasing shear flow involving a circular particle between two parallel 
plates with the 1024 x 256 mesh. The conditions are the same as those of figure 17. The colored filled circles are distributed 
to demonstrate the rotation. 



To directly demonstrate whether the reversibility can be captured, the distributions of the tracers for four consecu- 
tive time instants are shown in figure 19. As depicted in figure 19(a), the tracers are initially seeded on the concentric 
circles inside the solid to demonstrate the local displacements inside the solid. The bilinear interpolation to the tracer 
location is applied to identifying its velocity, and its position is temporally updated in a Lagrangian way. Figure 
19(b) shows the tracer distribution at the most deformed instant t = 4 when the particle is under the tank-treading 
like motion. After the wall velocities is set to be zero at t = 4, the tracer particles gradually move back toward the 
initial concentric circles with time. It should be noted that because the degree of freedom corresponding to the rigid 
rotation is allowed, the tracer distributions in figure 19(c) (d) turn in the clockwise directions about 80 degrees with 
respect to the initial distribution in figure 19(a). At the instant t = 8, when the same period as the shear-imposing 
stage (four unit time) has elapsed after the walls stop, the discrepancy between the tracer location and the concentric 
circle is clearly shown in figure 19(c), indicating that the recovery in the particle shape is still underway. After a 
sufficiently long time (t = 200), the tracers are found to be back in the concentric circles as shown in figure 19(d). We 
may say that the present Eulerian approach can capture the reversibility in shape under certain right circumstances, 
which will be discussed later. 

To assess the effect of the spatial resolution, the outlines of the fluid-structure interface, which are identified as the 
isolincs at <j) s = 1/2, for various number of grid points for different materials are shown in figure 20. As the spatial 
resolution is increased, the deformed shape more settles to a convergent curve at the instant t = 4 when the particle 
exhibits the most deformed shape, and the outline at t = 200 shows better recovery to the initial curve at t = 0. 

To further assess the grid convergence behavior in the particle deformation, the deformation mode is here inves- 
tigated. As explained in §IVB1, the deformation modes are determined using (54), (55) and (58). Due to the 
symmetry of the system, the odd-number-order modes R c ,2n+i and i? s ,2n+i arc identically zero. Temporal evolutions 
of the modal amplitudes \R n \ (n = 0, 2, 4 and 4) of the particle deformation for various number of grid points are 
shown in figure 21 (for the linear Mooney-Rivlin material) and in figure 22 (for the incompressible Saint Vcnant- 
Kirchhoff material). Both figure 21 and figure 22 demonstrate the convergence behavior of the profiles with increasing 
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FIG. 20: Outlines of the fluid-structure interface in the imposing-releasing shear flow involving a circular particle between two 
parallel plates for various number of grid points (N x x N v = 6-i x 16, 128 x 32, 256 x 64, 512 x 128 and 1024 x 256). The 
imposing-releasing shear scheme, the geometry, and the fluid properties are the same as those of figure 17. The left panels: the 
linear Mooney material with = 0, ci = 4, C2 = 2 and C3 = 0. The right panels: the incompressible Saint Venant-Kirchhoff 
material with fi 3 — 0, A£ ami = 6 and = 4. The upper panels: at t = 4. The lower panels: at t — 200. 



the number of grid points. After the walls stop at t = 4, the zeroth-order mode \Ro\ approaches 0.75, corresponding 
to the unstressed radius, in the case of the spatial resolution N x x N y = 256 x 64 or higher. However, the higher-order 
amplitudes |i?2|, \Ri\ and |i?6| at the fully developed stage obviously settle to some non-zero values. This tendency 
is more pronounced in the lower grid resolution, indicating that some spurious deformation remains. However, the 
deviation from zero for \R n \ (n ^ 0) vanishes exponentially as the spatial resolution is increased. 

To quantify the spurious residual deformation, the modal amplitudes \R n \ {n = 2, 4 and 6) at t = 200 as a function 
of the number of grid points are plotted in figure 23. The residual amplitudes are nearly proportional to iV" 1 , 
indicating the first-order accuracy in capturing unstressed shape. As demonstrated in §IVA for the much simpler 
system consisting of the fluid-structure layers, the present fluid-structure couping method involves the first-order 
accuracy, which is also reflected on the grid convergence of the reversibility in shape. 



2. Shape reversibility of a rectangular particle 

We also perform a reversibility test for a rectangular particle with a dimension of 2.375 x 1 to demonstrate the 
applicability of the method to an object with a larger aspect ratio and sharp corners. The initial setup is depicted 
in figure 24(a). Figure 24 visualizes the particle deformation and the velocity and vorticity fields for six consecutive 
time instants. Similar to the system involving the circular particle in figure 17, the elastic wave propagation and its 
attenuation are observed in figure 24(b)(c)(d). However, unlike the tank-treading motion in figure 17(e), the vorticity 
inside the particle at t = 4 is not entirely negative in figure 24(e), indicating that the particle does not experience 
the tank-treading or tumbling motion. It is because the rotational motion is geometrically suppressed due to the 
hydrodynamic interaction between the particle and the wall. As shown in figure 24(e), the left-top and right-bottom 
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FIG. 21: Time history of the n-th order modal amplitude \R n \ of the particle deformation in the imposing-releasing shear flow 
involving a circular particle between two parallel plates for various number of grid points (N x x N y — 64 x 16, 128 x 32, 256 x 64, 
512 x 128 and 1024 x 256). (a) n = 0, (b) n = 2, (c) n = 4, and (d) n = 6. The solid obeys a linear Mooney-Rivlin law with 
fi s = 0, ci = 4, C2 = 2 and C3 = 0. 



corners of the object are largely deformed. From the subsequent results on the shape reversibility in figure 26 and the 
grid convergence behavior in figure 27, we strongly envisage that large deformations arc resulted from the physical 
mechanism, not induced by numerical errors. After the shearing force is released at t = 4, the deformed particle 
gradually recovers the unstressed shape as shown in figure 24(f) (t = 6). 

Figure 25 shows the budget of the kinetic energy transport. Similar to the results in figure 18, the numerical error 
is much smaller than the variation of the contributions of the individual terms in (60), indicating that the energy 
exchange between the fluid and solid phases via the solid deformation is reasonably guaranteed. 

Figure 26 shows the tracer distributions for four consecutive time instants. After the particle well deforms at t = 4, 
the recovery of the material points toward the initial configuration is demonstrated. Even though the rotational 
motion of the particle is suppressed, the object at fully developed state in figure 26(d) slightly turns in the clockwise 
direction as compared with the initial distribution in figure 26(a). It would be remarkable to note that, though the 
left-top and right-bottom corners of the object are strongly stretched at t = 4, the object gradually resumes the shape 
of the original corners as the time goes on. 

Figure 27 shows the outlines of the fluid-structure interface for various number of grid points and different materials. 
The edges of the rectangle are obviously smeared out, which would be the inevitable effect of the numerical dissipation 
involved in the fifth-order WENO scheme, which is applied to updating <f) s . Nevertheless, similar to figure 20, with 
increasing the number of grid points, the particle shapes at t = 4 and t = 200 converge, and the reversibility in shape 
can be better attained. 
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FIG. 22: Same as figure 21, but the solid obeys an incompressible Saint Venant-Kirchhoff law with fi s = 0, XI 
Ml™* = 4. 
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V. CONCLUSION AND PERSPECTIVES 



A full Eulcrian simulation method for solving Fluid-Structure Interaction (FSI) problems has been developed. A 
volume-of-fluid formulation [23] was applied to describing the multi-component geometry. The temporal change in 
the solid deformation was described in the Eulcrian frame by updating a left Cauchy-Grecn deformation tensor, which 
was used to express the nonlinear Mooncy-Rivlin constitutive law. The validity of the present simulation method 
was established through comparisons with the analytical solution of the the oscillatory response in fluid-solid parallel 
layers, and also with the available simulation data of the solid motion in the lid-driven cavity flow [97] and the 
two-particle interaction in the Couette flow [12]. We confirmed that the present Eulerian approach can capture the 
reversibility in shape as long as the grid resolution is sufficiently high. Further, we demonstrated that the numerical 
accuracy due to the fluid-structure coupling is of the first-order with respect to the grid size. 

The significance of the present full Eulerian simulation method may be that the approach showed a feasibility of 
reducing the FSI coupling problem to a simple incompressible fluid flow solvers. Thus, the conventionally-used efficient 
computational techniques, such as the fast Fourier transform, and multi-grid method, are applicable. The present 
Eulcrian method is proved to be well-suited for using the voxel-based multi-component geometry on the fixed Cartesian 
system. Once the initial field of the solid volume fraction is given over the entire domain, the present Eulerian method 
enables one to carry out the FSI simulation without mesh generation procedure. The method promises to extend 
the possibility of the FSI simulation to certain additional classes of problems in the medical field, owing to a facility 
in incorporating the voxel data directly converted from medical images. The practical demonstration is the future 
subject of the present authors. 

To improve the accuracy in the present fluid-structure coupling to a level available for practical applications, it is 
important to capture the interface more sharply. We now use the fifth-order WENO method for advecting the solid 
volume fraction field, which temporally makes the interface numerically diffusive. As frequently used in the multiphase 
flow simulation, to suppress the numerical diffusion, elaborated techniques for the sharp interface advection such as 
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FIG. 23: Residual modal amplitudes \R n \ of the particle deformation at t = 200 versus the number N x of grid points in the 
imposing-releasing shear flow, (a): for the results in figure 21 employing the linear Mooney-Rivlin material (b): for the results 
in figure 22 employing the incompressible Saint Venant-Kirchhoff material. 



SLIC [52], PLIC [17, 89], and THINC [86] methods would be applicable. As an alternative of the VOF function, the 
level set function [56, 64] is another option. On the dynamic interaction, we now write the stress in a fluid-structure 
mixture form. Although the strain rate has a discontinuity across the fluid-structure interface, it is smoothed out 
at the grid scale in the present simulation method. The ideas of the immersed interface treatment [41, 42] and the 
localized strain formulation [53] would be effective to improve the accuracy in the fluid-structure coupling. It is a 
challenging task to overcome the multiphysics difficulty particularly associated with the difference in constitutive laws 
for fluid and solid. Improved accuracy in capturing the interface and robust time advancement [32, 33] are the ongoing 
subject of the present authors. 
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FIG. 26: Material point distribution in the imposing-releasing shear flow involving a rectangular particle between two parallel 
plates, with the 1024 x 256 mesh. The conditions are the same as those of figure 24. 
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FIG. 27: Outlines of the fluid-structure interface in the imposing-releasing shear flow involving a rectangular particle between 
two parallel plates for various number of grid points (N x x N y = 64 x 16, 128 x 32, 256 x 64, 512 x 128 and 1024 x 256). The 
imposing-releasing shear scheme, the geometry, and the fluid properties are the same as those of figure 24. The left panels: the 
linear Mooney material with /j, 3 — 0, C\ = 4, C2 = 2 and C3 = 0. The right panels: the incompressible Saint Venant-Kirchhoff 
material with fi s = 0, A£ am6 = 6 and /i£ am<i = 4. The upper panels: at t = 4. The lower panels: at t = 200. 
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Appendix A: Finite difference descriptions 



1. For mass conservation equation (11) and derivatives of the incremental pressure involved in (31) — (34) 



Using the operators in (38), we describe 

5i(v x )U,j , Sj(v v )\ 



(V.«) w = ^^ + ^™i. (Al) 



(3r¥>)i+J 



s j(<P)\ij+i 



(A2) 



5i{tp)\ t+ i, - 6i(tp)\ 



^2 



A 2 



(d x d y <p) i+ i tj+} 



<Pi+i,j+i ~ ~ <Pi+ij + <Pi,j 

A X Ay 



(A3) 



2. For momentum conservation equation (26) 



Here, we show the discretization for each term involved only in the ^-momentum equation. The permutations i <-> j 
and x O y lead to the corresponding discretization in the y-momentum equation. For a quantity q, we here introduce 
interpolation operators denoted by overlines such as 

t\i,j= +2 ' J 2 2 ' J , q 3 \i,j= J+2 2 2 - (A4) 

The advection terms [38]: 

v^\i,jSi(v x )\i t j +v x ' l \i + i, 3 Si(v x )\i + i^ 



(v x d x v x ) i+ i d 

v^ 3 \ i+ L tj ^5j{v x )\ i+ x tj _x + v^ J \ i+ i tj+ iS j (v x )\ i+ i !j+ i 



(vyd y v x )i. 



The pressure gradient and the divergence of the deviatoric stress tensors: 



(d x p) i+ ij 



L± x L\y 

where 
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(d- xy ) i+ i tj+ i = ([if + (fj, s - (J,f)(f> s \ i+ i !j+ i)(L xyii+ i d+ i + L yxi+ i + ((t> s a sh ,xy)i + i,j+i, 

{4>sPsh,xx)i,3 =|(2Cl - Ucs^l.j + ( 2C2 + 4c 3) tr (^)»J " - , '^/>'-.-..; I /!,•,.,., - 2c 2 -B^| 

(c/>sa sh , X y) i+hj+h =|(2c x + 2c 2 - 12c 3 )4| i+ i, i+ x + 4c 3 tr(B)| . +|| . + , js^+i^i. 

Considering S zz = 1, we write the trace of B as 

~ i 
tr(B)jj- = B X x,i,j + Byy i j + (t>s,i,j' 

3. For the advection terms in (13) and (22) 

For a quantity g (corresponding to -B X xj or -Byy) defined at the cell centroid («, j), we apply the fifth-order 
WENO scheme [36, 45] to the advection terms in (13) and (22). The advection term v x d x q is written as 



(v x d x q)i 
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where e is a positive tiny number to avoid division by zero, and 
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Likewise, (^yC^g)]^™ is computed using the interpolated advection velocity v y \ij. For 5^ defined at the cell apex 
(i + 5 > i + 2)1 usm S the interpolated velocities 1 j+i and v y \ i+ i j + ^, we evaluate v • VB xy in a similar manner. 



Appendix B: Spectral algorithm to And sharp interface solution for the parallel layers problem 

We here explain the sharp interface approach to solve the one-dimensional fluid-structure coupling problem by 
means of (pseudo) spectral method. We obtain accurate solutions, which are used for validating the present full 
Eulerian model by comparisons. 
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Due to the symmetry of the system with respect to y — illustrated in figure 5, we consider the upper half region 
y > and write the fluid and solid velocities Vf, v s in a Fourier series form 

v f (y,t) = Vj(t) + ±{V w (t) - Vj(t)) + «/,*(*) sin (Bl) 



(B2) 



fc=i 



where Vj is the velocity at the fluid-structure interface (y = L s ), Vw is the given velocity of the upper wall, Vf k and 
v s .k are expansion coefficients, and y = y — L s . The expressions (Bl) and (B2) satisfy the continuity of the velocity 
(vf — v s ) at the interface y = L s , the no-slip condition (vf = Vw) on the upper wall y = L a + Lf , and the symmetric 
condition (v s — 0) at y — 0. From (B2), we readily find the solid displacement u s as 



u 3 (y,t) 



Ui(t)y 
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where Ui and u Si k yield 
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From the momentum equations (2), with the stress expressions (10) and (17), we obtain 
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where <tnl denotes the nonlinear contribution in the solid stress with respect to the displacement. The definition of 
<7nl and the relation with the expansion coefficients ONL.fc are 



onl = 4c 3 



du s 
dy 
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From the orthogonality in the sine function, (B5) and (B6) are reduced to the modal relations 
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for 1 < k < 00. The continuity of the shear stress at the interface y = L s is 

lifiy w -Vi) 2(c 1 + c 2 )U I 
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The equation set to be solved consists of (B7) — (B10). In the numerical determination of the coefficients Vf u Si /-, Vj 
and E/j, we truncate the number of the modes appeared in the infinite series summation of (B10) up to k = K — 1. If 
K is chosen as an integer power of 2, the fast Fourier sine transform can be applied to efficiently evaluating Vf and u s 
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respectively given in (Bl) and (B3), and the fast Fourier cosine transforms determine the nonlinear part of the solid 
stress in a pseudo-spectral way 

8c 3 \ u i , nnUs - 1 ^Kj + |) 1 nk(J + |) 

CTNL ^^TTMg\i: + S^r C0S ^^| cos ^^' (B11) 

where 5 is the Kronecker delta. 

In the case of C3 = (the linear Mooney-Rivlin material), the system is linear since ctnl vanishes. Considering 
the wall velocity is Vyy(t) — Im(Vw exp(iwt)), we may apply the separation of variable to the velocities and the 
displacement 

Vi(t) =Im(Vf exp(iut)), 
v ft k(t) =Jm(vf t k exp(iwt)), (B12) 
«*,fc(t) =Im(ti S) fcexp(iwt)) ) 

which reduce the differential equations (B8) and (B9) with respect to £ into the algebraic ones. We readily find the 
expansion coefficients 
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In the case of C3 ^ (e.g., the incompressible Saint Venant-Kirchhoff material), the system is nonlinear since 
cnl 7^ 0, and thus the numerical time integration is needed. Here, it is carried out using the second-order Adams- 
Bashforth and Crank-Nicolson schemes. We here put an superscript (n) to a quantity to indicate the n-th time level 
(t = n(A£)). If all the variables at the n-th and (n — l)-th time levels are known, together with the prescribed wall 
velocity 

V$ +1) = V w ((n + l)(At)) = Im(VW exp(iu;(n + l)(At))), 
we update Ui, t>/,fc, u s ,k, and Vi at the (n + l)-th time level: 

= f** 2;2 , (B17) 
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where 
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After a sufficiently long computation, wc obtain temporally periodic solutions. 

We checked the convergence of the solution as a function of the truncated mode K. We confirmed that within the 
parameter range shown in figure 9, the results with K = 2048 arc accurate enough to be regarded as the reference 
solutions for comparison. 
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